Upregulation of IFNɣ-mediated chemokines dominate the immune transcriptome of muscle-invasive urothelial carcinoma

Tumor inflammation is prognostically significant in high-grade muscle-invasive bladder cancer (MIBC). However, the underlying mechanisms remain elusive. To identify inflammation-associated immune gene expression patterns, we performed transcriptomic profiling of 40 MIBC archival tumors using the NanoString nCounter Human v.1.1 PanCancer Panel. Findings were validated using the TCGA MIBC dataset. Unsupervised and supervised clustering identified a distinctive immune-related gene expression profile for inflammation, characterized by significant upregulation of 149 genes, particularly chemokines, a subset of which also had potential prognostic utility. Some of the most enriched biological processes were lymphocyte activation and proliferation, leukocyte adhesion and migration, antigen processing and presentation and cellular response to IFN-γ. Upregulation of numerous IFN-γ-inducible chemokines, class II MHC molecules and immune checkpoint genes was detected as part of the complex immune response to MIBC. Further, B-cell markers linked to tertiary lymphoid structures were upregulated, which in turn is predictive of tumor response to immunotherapy and favorable outcome. Our findings of both an overall activated immune profıle and immunosuppressive microenvironment provide novel insights into the complex immune milieu of MIBC with inflammation and supports its clinical significance for predicting prognosis and immunotherapeutic responsiveness, which warrants further investigation. This may open novel opportunities to identify mechanisms for developing new immunotherapeutic strategies.

www.nature.com/scientificreports/ and the underlying molecular pathways are still lacking. Future progress in MIBC therapeutic approaches will require a strong understanding of the immunology of this disease. Therefore, research into the immune landscape of MIBC may lead to the development of more precise prognostication, patient stratification and personalized treatment decisions. In this regard, transcriptomic data is increasingly being used to determine the tumor immune microenvironment and these methods have been shown to be very efficient and reliable 17 . Targeted gene expression panels, such as NanoString panels, have the potential to assess the tumor microenvironment and immune milieu of MIBCs comprehensively and ascertain the inflammatory status of a tumor by quantifying chemokines, cytokines, cell surface proteins and transcription factors along with immune checkpoints.
In the current study, we investigated gene expression profiles of immune-related genes (IRGs) as detected by NanoString technology to generate knowledge about the tumor immune milieu in highly inflamed versus low-inflamed MIBCs and better understand the potential immunobiological mechanisms for the association of inflammation with clinical outcome.

Results
Sunnybrook high-grade urothelial carcinoma cohort. The clinicopathologic characteristics of the MIBC cohort are listed in Table 1. Patients with high inflammation intra-tumorally also had high inflammation at the invasive front. There was a significant association of tertiary lymphoid structures (TLS) with high inflammation (19/20) and only one case of low inflammation containing TLS (P < 0.0001, Fig. 1). Four low inflammation tumors had small lymphoid aggregates lacking defined germinal centers. These were remote from the tumor edge in the perivesical adipose tissue.
Global immune gene expression patterns of high-grade urothelial carcinomas. To characterize inflammation-associated immune gene expression patterns in high versus low inflammation MIBC, NanoString profiling coupled to the nCounter human PanCancer immune panel was utilized to interrogate the mRNA Table 1. Clinicopathologic characteristics of the Sunnybrook NanoString cohort. CIS Carcinoma In Situ, pT Pathologic T stage, N/A not available.  www.nature.com/scientificreports/ expression profiles. In total, expression of 730 mRNAs from 24 different immune cell types related to both innate and adaptive immune responses were screened. Of these, 532 (72.9%) IRGs were expressed above background threshold levels. Unsupervised hierarchical clustering of gene expression levels using Euclidean correlation segregated the 40 MIBCs into two clusters, representing different immune profiles ( Fig. 2A-B). Further, K-Means consensus clustering analysis also yielded the same two clusters, suggesting that two robust immune expression subtypes exist in our cohort. Cluster I consisted of primarily high inflammation cases (15/17, 88.2%) while cluster II consisted primarily of cases with low inflammation (18/23, 78.3%).
Differentially expressed immune related genes associated with inflammation. To investigate the inflammation-associated modulation in immune gene expression patterns between high and low inflammation, linear regression analysis with multiple-testing correction was utilized. As expected, the immune signature of high inflammation tumors was characterized by significant (FDR P ≤ 0.05) over twofold altered expression of 149 genes associated with key immune response pathways, 99.3% of which were upregulated. The top 25 most differentially expressed genes are listed in Table 2 and these included the TLS marker CXCL13, B-cell related genes (CD19, CD79, TNFRSF17, POU2AF1 and PIK3CD), T-cell lineage and T-effector genes (CD2, CD3D, CD3E, CD6, CD4, CD8, IL2RG and CD44) as well as genes expressed on T-regulatory cells (FOXP3, CTLA4 and ICOS). The leukocyte common antigen, CD45, was also significantly upregulated in high inflammation tumors. Assessment of the expression levels of B-cell, T-cell, CD8 T-cell, Th1 cell and dendritic cell markers also showed upregulation of these genes in high inflammation tumors (Fig. 3A).
Lastly, the immune checkpoint genes (ICGs) CTLA-4, PD-L1, PD-2, IDO1 and LAG3 were significantly upregulated in high inflammation tumors while PD1 expression was below threshold of detection in tumors with low inflammation.

Pathway analysis.
To determine the inflammation-associated immunobiological pathways governing the differences between low and high inflammation, g:Profiler functional enrichment analysis was utilized. At FDR www.nature.com/scientificreports/ corrected P ≤ 0.05, we found that upregulated transcripts were significantly enriched for 85 GO biological processes, 40 KEGG pathways and 21 Reactome pathways (Supplemental Fig. 1). The most enriched biological processes distinguishing high from low inflammation were related to lymphocyte activation and proliferation, leukocyte adhesion and migration, cellular response to tumor necrosis factor, antigen processing and presentation, cellular response to IFN-γ, neutrophil activation, positive regulation of ERK1 and ERK2 cascade and calcium ion transport. The top overrepresented pathways included cytokine-cytokine receptor interaction, hematopoietic cell lineage, signaling by Interleukins, Th1, Th2 and Th17 cell differentiation, chemokine signaling, neutrophil degranulation, costimulation by the CD28 family, NF-kappa B signaling pathway, Jak-STAT signaling pathway, antigen processing and presentation and T cell receptor signaling pathway (Supplemental Fig. 1).
Prognostic Immune-related mRNAs. Cox proportional hazards model was used to assess whether DEG between tumors with high and low invasive front (FDR-corrected P ≤ 0.001, n = 71) are prognostically significant for EFS. To test this effect, patients were median-dichotomized to low and high expression based on the mRNA abundance levels of each DEG. The list of genes, hazard ratios, 95% confidence intervals (CI), P values and www.nature.com/scientificreports/  www.nature.com/scientificreports/ BH corrected p-values (adjusted for multiple testing) are listed in Table 3. We identified 13 genes differentially expressed between high and low inflammation tumors as significant potential predictors of EFS.
TCGA data. To validate the gene expression signatures we identified in association with high and low inflammation, we performed an immune transcriptomic analysis of the immune gene expression profiles of 730 genes derived from the NanoString™ panel utilized for our Sunnybrook cohort. The TCGA cohort consisted of 411 chemotherapy-naïve, muscle-invasive, high-grade urothelial tumors. Unsupervised hierarchical clustering of gene expression levels using Euclidean correlation segregated the 411 TCGA high-grade urothelial carcinomas into two large clusters, representing different immune profiles ( Fig. 2C-D). Cluster I could also be subdivided into two sub-clusters (Cluster Ia and Ib). KMeans consensus clustering analysis revealed that the TCGA cohort can be more robustly segregated into two compared to three clusters (Fig. 2C). Cluster I consisted of approximately equal distribution of cases with and without inflammatory infiltrate (146/287 49% and 141/287 51%, respectively). This cluster can be further subdivided into Cluster Ia, consisting primarily of cases without inflammatory infiltrate (74/102, 73%) and Cluster Ib; 39% (72/185) without and 61% (113/185) with inflammatory infiltrate. Cluster II consisted primarily of cases with inflammatory infiltrate (106/124, 85%).
Linear regression analysis with multiple-testing correction revealed that cases with inflammatory infiltrate were characterized by significant (FDR P ≤ 0.05) over twofold altered expression of 221 genes associated with key immune response pathways, 98.6% of which were upregulated. Among these, 112 were also significantly (FDR P ≤ 0.05) over twofold altered in our NanoString Sunnybrook cohort. Additional 32 genes were significantly (FDR P ≤ 0.05) altered in both Sunnybrook and TCGA cohorts, but did not reach twofold change in the latter. Therefore, 96.6% (144/149) of DEGs in our NanoString Sunnybrook cohort were validated in the TCGA cohort. Comparison of the fold change and FDR P-values for the top 25 mRNAs is presented in Table 2.
Functional enrichment analysis with g:Profiler validated that immunobiological pathways governing the differences between cases with and without inflammatory infiltrate in the TCGA cohort are largely similar to the pathways that differentiate between low and high inflammation in the Sunnybrook NanoString cohort (Supplemental Fig. 1).

Discussion
Multigene immune signatures represent a robust means of capturing the key players in MIBC tumor-associated inflammation. In the current study, we used a NanoString immune profiling panel to analyze and characterize inflammation-associated immune gene expression patterns in high inflammmation versus low-inflammation MIBCs. Unsupervised hierarchical clustering analysis identified two clusters characterized by distinctive immune profiles. We found that those of the high-inflammation cluster had a higher expression of most immune cell markers, including Th1 cells, CD8 + and B cells. Further, numerous chemokines, which play a crucial role in modulating anti-tumor immune response, were also up-regulated in high inflammation tumors. Conversely, the low inflammation cluster had a relatively lower abundance of immune cell markers. We validated our findings and obtained similar results analyzing the TCGA MIBC dataset. Analyzing DEG between low-and high-inflammation tumors in our cohort, we identified a signature of 149 inflammation-associated IRGs, a subset of which also had a good ability to predict EFS.
Upregulation of CCL4, CCL5, CXCL9, CXCL10 and CXCL11 in high inflammation MIBC is noteworthy, as their presence in tumors promotes T cell infiltration, a key mechanism in antitumor immunity and could explain improved patient outcome 18 . The upregulation of these chemokines is induced by IFN-γ, secreted by Th1 cells 19 . Indeed, Th1-associated markers were upregulated in high inflammation tumors. Moreover, expression Table 3. Univariate Cox regression analysis of differentially expressed genes between with high-versus lowinvasive front for predicting EFS in muscle invasive urothelial carcinomas treated by cystectomy. HR Hazard Ratio, CI Confidence Interval, BH-FDR Benjamini-Hochberg procedure to decreases the false discovery rate. www.nature.com/scientificreports/ of a number of cytotoxic markers induced by IFN-γ, namely granzymes A, B and K as well as PRF1 was also elevated in high inflammation MIBC, consistent with the presence of activated cytotoxic T cells together with Th1 cells in the tumors of these patients 9 . Further, we detected the upregulation of numerous additional IFNγ-inducible chemokines and class II MHC molecules associated with anti-tumor inflammation and improved patient outcome 20 . This is in agreement with previous studies that have shown that "hot" tumor microenvironment characterized by T cell infiltration and pro-inflammatory cytokine production present higher response rates to immunotherapy compared to infiltration-excluded or "cold" tumors 21 . Further, it has been reported that patients with high immune infiltration and IFN-γ-related mRNA expression signature correlated with a clinical benefit 22 . Of note, the IFN-γ-inducible chemokines CXCL9, CXCL10 and CXCL11 activate the chemokine receptor CXCR3, a variant transcript of which has been associated with response to neoadjuvant chemotherapy 23 . It is proposed as one of the potential mechanisms underlying the increased responsiveness of the immune-infiltrated basal subtype of MIBC to neoadjuvant chemotherapy 24 . However, the expression of CXCR3 was not linked to inflammation in our cohort, possibly because we did not analyze the expression of the individual CXCR3 isoform subpopulations 25 . Further research into the relevance of CXCR3 isoforms relative to the inflammatory tumor milieu for the response to chemotherapy is needed to elucidate the role of the immune system in the response to chemotherapy in MIBC.
High inflammation MIBC tumors also showed significant upregulation of the IFN-γ-inducible ICGs IDO1 and PD-L1 as well as the ICGs CTLA-4 and LAG3. The stimulation of these ICGs likely facilitates MIBC protumor immune tolerance, suggesting that immune checkpoint blockade therapies could be more suitable for MIBC patients with high-compared to low-inflammation. This is currently being investigated in a phase II clinical trial titled DUTRENEO of durvalumab and tremelimumab versus chemotherapy as a neoadjuvant approach to MIBC treatment of patients selected based on IFN-γ immune signature 26 .
In addition to IFN-γ inducible genes, our analysis identified that CXCL13, a B-cell chemoattractant, and CD79A, part of the B-cell receptor signaling complex, are among the top most differentially expressed genes between high-and low-inflammation MIBCs. Expression of both CXCL13 and CD79A has been linked to TLS, which in turn is predictive of tumor response to immunotherapy and favorable outcome in many tumor types, including MIBC [27][28][29][30] . Importantly, recent data reported a correlation between CXCL13 expression and response to ICIs 27,28 . In agreement with these previous studies, there was a significant association of CXCL13 and TLS with high inflammation in our cohort. Further, recent studies have shown that TLS promotes recruitment of lymphocytes to the tumor microenvironment by expressing chemokines such as CXCL10, CCL19 and CCL21 31 . Consistent with previous findings, these chemokines were all upregulated in high inflammation MIBCs in our cohort. Taken together, this could explain improved outcome in patients with highly inflamed phenotype.
TLS has been closely associated with the recruitment, activation and proliferation of B cells 32 . Consistent with this finding, our analysis has identified upregulation of numerous B cell markers in MIBC with high inflammation. Depending on the B-cell subtype, mode of activation and the microenvironment, these cells can either contribute to upregulation of T-cell responses or they can exert immunoregulatory functions and participate in the downregulation of T-cell immunity. Recently, CD19 + tumor-infiltrating B-cells were identified as an independent favorable prognostic factor in MIBC and serve as antigen-presenting cells to activate CD4 + T cell in the tumor microenvironment 33 . In agreement with these findings, we found upregulation of CD19 and CD4 in high inflammation MIBC. Further, antigen processing and presentation was one of the most enriched biological processes distinguishing high from low invasive front inflammation in our study.
The current findings add to a growing body of literature characterizing the immune genes and signaling pathways in MIBC utilizing the TCGA MIBC cohort. In 2017, Ren R. et al. 34 performed transcriptomic analysis of immune genes in the TCGA MIBC cohort in association with MIBC molecular subtypes (TCGA subtypes 2014) 34 . They found an increased expression of immune-associated genes in Cluster IV and underactive immune environment in Cluster I, which could be of potential significance for immunotherapies 34 . The study also identified responses to IFN-γ as one of the most common overrepresented pathways in the immune transcriptome of MIBCs and investigated it in the context of TCGA molecular subtypes. Zhou et al. explored 29 immune gene signatures in MIBC TCGA data to identify 3 MIBC immune subtypes characterized by high, medium and low immune gene expression signatures 35 . The group then used a prediction algorithm to estimate the level of immune cell infiltration, stromal content, and tumor purity of each MIBC. They found that immunity high subtype had a high level of immune infiltration, proportion of specific tumor infiltrating lymphocytes, HLA abundance, PD-L1 expression levels, were more likely to respond to ICIs and had relatively better clinical outcomes 35 . Lastly, a pan-cancer study based on TCGA data identified 6 immune subtypes for 33 non-hematologic tumors, including MIBC 36 . They also used prediction algorithms to estimate immune cell content from the TCGA multiomics data among the 6 immune subtypes to identify immunogenomic features that were predictive of outcome. This study found the immunosuppressed subtype had the worst prognosis while IFN-γ dominant and inflammatory subtypes had the most favorable prognosis. Both studies comprehensively studied the immunogenomics characteristics of MIBC TCGA cohort and suggested a prognostic value for immene signature based molecular subtyping of MIBC. This is distinct from our current study, which characterized the immune transcriptome signatures associated with the degree of inflammation in MIBC as seen on histopathological review.
In this study, we used retrospective FFPE material from a small cohort and the TCGA dataset, which has been generated from frozen bladder tumors. Another possible difference between the cohorts is in the morphological features used to define inflammation. Within the Sunnybrook cohort, five patients received remote intravesical Bacille Calmette Guerin (BCG). It is unclear what (if any) impact this had on the overall immune composition. Despite these limitations, the IRG transcriptional landscape was similar across the two cohorts. The findings reported in the Sunnybrook cohort are based on a limited sample size. Therefore, further research with larger sample sizes are needed to validate the immune signature and more experimental research is needed www.nature.com/scientificreports/ to fully elucidate the biological mechanisms underlying inflammation, response to immunotherapy, response to chemotherapy and clinical outcome in MIBC. We discerned two distinct transcriptomic immune profiles for high-and low-tumor inflammation MIBCs and identified a 149 IRG signature associated with high-tumor inflammation, a subset of which also have potential prognostic utility. The immune profile of high invasive front inflammation was characterized by both upregulation of pro-inflammatory chemokines and immunosuppressive checkpoints. Notably, we found upregulation of IFNɣ-mediated chemokine signature in high inflammation MIBC as well as elevated PD-L1 expression. These findings were validated in the TCGA dataset. The findings further enhances our understanding of the tumor immune milieu and supports its clinical importance for predicting prognosis and immunotherapeutic responsiveness in MIBC, which warrants further investigation.

Materials and methods
Sunnybrook cohort clinicopathological variables. Institutional approval was obtained from the Sunnybrook Health Sciences Centre research ethics board for the study (REB 187-2016) and the need for written informed consent for this retrospective study was waived under the approval of the IRB. All procedures were performed in accordance with the Declaration of Helsinki and its later amendments or comparable ethical standards.
Retrospective formalin-fixed paraffin-embedded (FFPE) cystectomy specimens from 40 high-grade muscleinvasive urothelial carcinomas treated at Sunnybrook Health Sciences Centre were included in the current study, 20 with high inflammation and 20 with low inflammation (see definition below). This is a subset of randomly selected cases from a larger cohort of 235 cases characterized in our previous study 6 . Clinicopathological and follow-up data ( Table 1) was obtained according to protocols approved by the Research Ethics Board of Sunnybrook Health Sciences Centre. All tumors had predominant urothelial histology. Event-free survival (EFS) was defined as the time from surgery to the following defined events: recurrence, death or palliative care. Subjects who did not experience an event were censored from date of last follow-up. MIBC patients in this cohort did not receive adjuvant or neoadjuvant immune checkpoint inhibitors. Two patients with high and three with low inflammation received neoadjuvant gemcitabine and cisplatin. Five patients received prior intravesical BCG: two in the high invasive front group and three in the low invasive front cohort. All H&E slides were reviewed by a urological pathologist (MRD) for the following: histology, grade (2004 WHO/International Society of Urological Pathology), stage (AJCC/TNM), LVI, margin status (a positive ST margin was defined as tumor extending to ink), carcinoma in situ (CIS) and lymph node metastases.

Assessment of inflammation in the Sunnybrook cohort.
For each case, slides containing urothelial carcinoma were selected and inflammation was semiquantitatively graded by two pathologists (AH/MRD) using the four-point system developed by Klintrup-Makinen 37 as previously described 6 . Inflammation was assessed at both the invasive front (defined as the deepest interface of carcinoma with stroma) and centrally within the tumor and assigned a score of 0, 1, 2 or 3 at each site. All inflammatory cells were included in the assessment. At the invasive front, inflammation was scored as follows: 0, none; 1, patchy infiltrate; 2, bandlike infiltrate; and 3, prominent with destruction of urothelial carcinoma cells. Only the inflammatory component in contact with and directly adjacent to urothelial carcinoma was scored. Within the tumor, inflammation was scored as follows: 0, absent; 1, scattered inflammatory cells; 2, small aggregates of inflammatory cells; and 3, dense intratumoral inflammation. The scores were averaged across both components and subsequently dichotomised as low (score 0 or 1) and high (score 2 or 3). The presence/absence of tertiary lymphoid structures (TLS) were assessed using published criteria 30 . We defined TLS as well defined clusters of lymphocytes including those with germinal centres. Loose lymphoid collections were recorded as lymphoid aggregates.
Sunnybrook cohort RNA extraction and NanoString Gene expression profiling. RNA was extracted from macrodissected FFPE tissue Sections (4 to 10 sections per case, each 5 µm thick) using the High-Pure FFPET RNA Isolation Kit (Roche Diagnostics, Indianapolis, IN, USA), according to the manufacturer's protocol. RNA was quantified with the Thermo Scientific Nanodrop 1000 (Thermo Fisher Scientific, Wilmington, DE). RNA profiling was performed with 250 ng (quantified using NanoDrop-1000, Thermo Fisher Scientific, Waltham, Massachusetts, USA) of RNA using the NanoString nCounter Human V.1.1 Pan-Cancer Immune Profiling Panel (NanoString Technologies Inc, Seattle, Washington, USA), according to the manufacturer's instructions.
Raw gene expression data was analyzed using NanoString's software nSolver V.4.0 with the Advanced Analysis 2.0 plugin. Data normalization was performed using internal negative control probes, synthetic positive controls and 36 selected housekeeping genes that were identified using the nSolver normalization module, which uses the geNorm algorithm (https:// genorm. cmgg. be/). Data was normalized using geometric mean, log2-transformed and then used as input for further analysis. KMeans consensus clustering was performed using Pearson distance and average linkage in GenePattern (Broad Institute of MIT and Harvard, http:// www. broad. mit. edu/ cancer/ softw are/ genep attern). For each value of k from 2 to 10, 1000 iterations were performed. The consensus cumulative distribution function (cdf) and Δ area plots were examined and k = 2 was determined to be the ideal model. Unsupervised hierarchical clustering analysis using Euclidean correlation was performed using the R package ComplexHeatmap version 1.12.0. P values were adjusted using the Benjamini-Hochberg (BH) method to control the false discovery rate (FDR). Statistically significant, differentially expressed genes (DEGs) were defined as those with expression levels corresponding to a log2 ratio > 1 or < − 1 and BH q < 0.05. Functional and pathway enrichment analysis was performed using g:Profiler (V.e94_eg41_p11_9f195a1) with a Benjamini-Hochberg (BH-FDR) multiple testing correction method applying significance threshold of 0.05 38  www.nature.com/scientificreports/ TCGA cohort. Level three RNA sequencing (RNA-seq) data of The Cancer Genome Atlas (TCGA) Urothelial Bladder Carcinoma and corresponding clinical information were downloaded from The Genomic Data Commons Data Portal (https:// portal. gdc. cancer. gov/) and log2-transformed prior to analysis. Accession numbers and HGNC IDs were obtained for each ENSEMBL Gene ID in the TCGA database using the BioMart online tool (www. ensem bl. org/ bioma rt). The 730 genes represented on the NanoString's Pan-Cancer Immune Profiling Panel were then selected for downstream analysis (by matching Accession numbers and/or HGNC IDs). Data regarding immune infiltrates was obtained from publication by Robertson et al. 39 Table S2: Supplemental data was downloaded from Robertson et al. 39 , data in tab S2.4 Immune infiltrate 39 was dichotomized (present/absent) and matched to RNA-seq data using TCGA_IDs.
Statistical analysis. All statistical analyses were performed with SPSS 24.0 (IBM Corporation, Armonk, New York, USA). Univariate Cox proportional hazards analysis was used to correlate gene expression with EFS. BH corrected two-sided P-values less than 0.05 were considered statistically significant.

Data availability
Data are available from the corresponding author on reasonable request.